C
C     ..................................................................
C
C	 SUBROUTINE TETRA
C
C	 PURPOSE
C	    COMPUTE A TETRACHORIC CORRELATION COEFFICIENT BETWEEN TWO
C	    VARIABLES WHERE DATA IN BOTH VARIABLES HAVE BEEN REDUCED
C	    ARTIFICIALLY TO TWO CATEGORIES.
C
C	 USAGE
C	    CALL TETRA (N,U,V,HU,HV,R,RS,IE)
C
C	 DESCRIPTION OF PARAMETERS
C	    N  - NUMBER OF OBSERVATIONS
C	    U  - INPUT VECTOR OF LENGTH N CONTAINING THE FIRST VARIABLE
C		 REDUCED TO TWO CATEGORIES
C	    V  - INPUT VECTOR OF LENGTH N CONTAINING THE SECOND VARIABLE
C		 REDUCED TO TWO CATEGORIES
C	    HU - INPUT NUMERICAL CODE INDICATING THE HIGHER CATEGORY OF
C		 THE FIRST VARIABLE.  IF ANY VALUE OF VARIABLE U IS
C		 EQUAL TO OR GREATER THAN HU, IT WILL BE CLASSIFIED AS
C		 THE HIGHER CATEGORY, OTHERWISE AS THE LOWER CATEGORY.
C	    HV - SAME AS HU EXCEPT THAT HV IS FOR THE SECOND VARIABLE.
C	    R  - TETRACHORIC CORRELATION COMPUTED
C	    RS - STANDARD ERROR OF TETRACHORIC CORRELATION COMPUTED
C	    IE - ERROR CODE
C		 0 - NO ERROR
C		 1 - UNABLE TO COMPUTE A TETRACHORIC CORRELATION DUE TO
C		     THE FACT THAT AT LEAST ONE CELL SHOWS ZERO FRE-
C		     QUENCY IN THE 2X2 CONTINGENCY TABLE CONSTRUCTED
C		     FROM INPUT DATA.  IN THIS CASE, R AND RS ARE SET
C		     TO 10**75.  (SEE GUILFORD, 1956)
C		 2 - THE ROOT SOLVER GIVES MULTIPLE ROOTS, OR NO ROOTS,
C		     R, IN THE INTERVAL (-1,1) INCLUSIVE. R AND RS ARE
C		     SET TO 10**75.
C		 3 - UNABLE TO COMPUTE A SATISFACTORY VALUE OF TETRA-
C		     CHORIC CORRELATION USING NEWTON-RAPHSON METHOD OF
C		     APPROXIMATION TO THE ROOT OF THE EQUATION.  R AND
C		     RS ARE SET TO 10**75.  SEE SUBROUTINE POLRT ERROR
C		     INDICATORS.
C		 4 - HIGH ORDER COEFFICIENT OF THE POLYNOMIAL IS ZERO.
C		     SEE SUBROUTINE POLRT ERROR INDICATORS.
C
C	 REMARKS
C	    VALUES OF VARIABLES U AND V MUST BE NUMERICAL, AND
C	    ALPHABETIC AND SPECIAL CHARACTERS MUST NOT BE USED.
C	    FOR A DEPENDABLE RESULT FOR TETRACHORIC CORRELATION,
C	    IT IS RECOMMENDED THAT N BE AT LEAST 200 OR GREATER.
C
C	 SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C	    NDTRI
C	    POLRT--THIS POLYNOMIAL ROOT ROUTINE WAS SELECTED BECAUSE OF
C		   ITS SMALL STORAGE REQUIREMENT.  OTHER SSP ROUTINES
C		   WHICH COULD REPLACE POLRT ARE PRQD AND PRBM.  THEIR
C		   USE WOULD REQUIRE MODIFICATION OF TETRA.
C
C	 METHOD
C	    REFER TO J. P. GUILFORD, 'FUNDAMENTAL STATISTICS IN PSYCHO-
C	    LOGY AND EDUCATION', MCGRAW-HILL, NEW YORK, 1956, CHAPTER 13
C	    AND W. P. ELDERTON, 'FREQUENCY CURVES AND CORRELATION' 4-TH
C	    ED., CAMBRIDGE UNIVERSITY PRESS, 1953, CHAPTER 9.
C
C     ..................................................................
C
      SUBROUTINE TETRA (N,U,V,HU,HV,R,RS,IE)
C
      DIMENSION XCOF(8),COF(8),ROOTR(7),ROOTI(7)
      DIMENSION U(1),V(1)
      DOUBLE PRECISION X31,X32,X312,X322
C
C	 CONSTRUCT A 2X2 CONTINGENCY TABLE
C
      A=0.0
      B=0.0
      C=0.0
      D=0.0
      DO 40 I=1,N
      IF(U(I)-HU) 10, 25, 25
   10 IF(V(I)-HV) 15, 20, 20
   15 D=D+1.0
      GO TO 40
   20 B=B+1.0
      GO TO 40
   25 IF(V(I)-HV) 30, 35, 35
   30 C=C+1.0
      GO TO 40
   35 A=A+1.0
   40 CONTINUE
C
C	 TEST WHETHER ANY CELL IN THE CONTINGENCY TABLE IS ZERO.
C	 IF SO, RETURN TO THE CALLING ROUTINE WITH R=0.0 AND IE=1.
C
      IE=0
      IF(A) 60, 60, 45
   45 IF(B) 60, 60, 50
   50 IF(C) 60, 60, 55
   55 IF(D) 60, 60, 70
   60 IE=1
      GO TO 86
C
C	 COMPUTE P1, Q1, P2, AND Q2
C
   70 FN=N
      P1=(A+C)/FN
      Q1=(B+D)/FN
      P2=(A+B)/FN
      Q2=(C+D)/FN
C
C	 FIND THE STANDARD NORMAL DEVIATES AT Q1 AND Q2, AND THE
C	 ORDINATES AT THOSE POINTS
C
      CALL NDTRI (Q1,X1,Y1,ER)
      CALL NDTRI (Q2,X2,Y2,ER)
C
C	 COMPUTE THE TETRACHORIC CORRELATION COEFFICIENT
C
      IF(X1) 76, 72, 76
   72 IF(X2) 76, 74, 76
   74 R=0.0
      GO TO 90
   76 XCOF(1)=-((A*D-B*C)/(Y1*Y2*FN*FN))
      XCOF(2)=1.0
      XCOF(3)=X1*X2/2.0
      XCOF(4)=(X1*X1-1.0)*(X2*X2-1.0)/6.0
      X31=DBLE(X1)
      X32=DBLE(X2)
      X312=X31**2
      X322=X32**2
      XCOF(5)=SNGL(X31*(X312-3.0D0)*X32*(X322-3.0D0)/24.0D0)
      XCOF(6)=SNGL((X312*(X312-6.0D0)+3.0D0)*(X322*(X322-6.0D0)+3.0D0)
     1	      /120.0D0)
      XCOF(7)=SNGL(X31*(X312*(X312-10.0D0)+15.0D0)*X32*(X322*(X322-10.0
     1	      D0)+15.0D0)/720.0D0)
      XCOF(8)=SNGL((((X312-15.0D0)*X312+45.0D0)*X312-15.0D0)*(((X322-
     1	      15.0D0)*X322+45.0D0)*X322-15.0D0)/5040.0D0)
C
      CALL POLRT (XCOF,COF,7,ROOTR,ROOTI,IER)
C
      J=0
      IF(IER) 78, 78, 84
   78 DO 82 I=1,7
      IF(ABS(ROOTI(I))-.5*ABS(ROOTR(I))*1.0E-6)79,79,82
   79 R=ROOTR(I)
      IF(ABS(R)-1.0)81,81,80
   80 R=1.E38
      GO TO 82
   81 J=J+1
   82 CONTINUE
      IF(J-1)83,88,83
   83 IE=2
      GO TO 86
C
C	 UNABLE TO COMPUTE R
C
   84 IE=IER
   86 R=1.0E38
      RS=R
      GO TO 100
   88 IF(R-1.0E38)90,83,83
C
C	 STANDARD ERROR OF R=0.0
C
   90 RS= SQRT(P1*P2*Q1*Q2)/(Y1*Y2* SQRT(FN))
C
  100 RETURN
      END

COF,COF,7,ROOTR,ROOTI,IER)
C
      J=0
